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Abstract: We present the results of our computation of the chiral condensate with 
Nf = 2 and Nj = 2 + 1 + 1 flavours of maximally twisted mass fermions. The condensate 
is determined from the Dirac operator spectrum, applying the spectral projector method 
proposed by Giusti and Liischer. We use 3 lattice spacings and several quark masses at 
each lattice spacing to reliably perform the chiral and continuum extrapolations. We study 
the effect of the dynamical strange and charm quarks by comparing our results for Nf = 2 
and Nf = 2 + 1 + 1 dynamical flavours. 
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1 Introduction 

One of the most important phenomena of QCD is the spontaneous breaking of chiral sym- 
metry. This purely non-perturbative phenomenon was subject to many analyses in Lattice 
QCD, using different fermion discretizations and methods. In particular, the chiral con- 
densate - the order parameter of spontaneous chiral symmetry breaking - can be extracted 
from chiral perturbation theory fits of the quark-mass dependence of light pseudoscalar me- 
son observables [1-9], the topological susceptibility [8, 10-13] or the pion electromagnetic 
form factor [14]. Other methods include using the e-regime expansion and/or chiral ran- 
dom matrix theory [7, 13, 15-22], Wilson chiral perturbation theory fits of the integrated 
spectral density [23-25] and a calculation directly from the quark propagator [26, 27]. A 
summary of recent results is provided in Ref. [28]. 

The chiral condensate is related to the spectral density of the Dirac operator via the 
Banks-Casher relation [29]: 

lim lim lim p{X,m) = — , (1-1) 
A^O m-i-O y-s-oo vr 

where A is the modulus of the eigenvalue, m the quark mass, p(A, m) the spectral density, 
V the volume and S is the chiral condensate in the infinite- volume and in the chiral limit. 
Clearly, the triple limit on the left-hand side of the above equation makes it impractical to 
evaluate the chiral condensate on the lattice. 

However, recently a method has been proposed [30] to effectively make use of the Banks- 
Casher relation and explore the chiral properties of QCD on the lattice, in particular to 
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compute the chiral condensate. The method has also other apphcations, e.g. it ahows to 
compute the topological susceptibility or renormalization constants. Briefly, the method 
consists in stochastically evaluating the mode number, i.e. the number of eigenmodes of 
the Dirac operator below some spectral threshold value and using the dependence of this 
number of eigenmodes on the threshold value to calculate the observable of interest. In the 
following, we will refer to this method as spectral projectors. One of its essential advantages 
for computing the mode number is the fact that it is very effective in terms of computational 
cost - the required computational effort grows linearly with the lattice volume instead of 
quadratically, as is the case for a direct computation of eigenmodes and counting their 
number below the spectral threshold. 

In this paper, we report our results for the chiral condensate with Nj = 2 and Nf = 
2 + 1 + 1 Wilson twisted mass fermions at maximal twist. Preliminary results of our 
computations for the Nf = 2 + 1 + 1 case were presented in Ref. [31]. The paper is 
organized as follows. In the second section we provide a short description of the spectral 
projector method. In section 3, we describe our lattice setup. Section 4 presents our results 
for the chiral condensate both in the A'^^ = 2 and the Nf = 2 + 1 + 1 case. In section 5, 
we summarize and compare with other determinations of the chiral condensate found in 
literature. An appendix presents our tests of the method. 

2 Spectral projectors and chiral condensate 

Many interesting properties of the chiral regime of QCD can be understood from the be- 
haviour of quantities related to the low-lying spectrum of the Dirac operator. One of such 
spectral quantities, essential in the determination of the chiral condensate, is the mode 
number, i.e. the number of eigenvectors of the massive Hermitian Dirac operator D^D with 
eigenvalue magnitude smaller than some threshold value M^. We will denote this mode 
number by u{M,^), where ^ is the quark mass. 

Here we provide a short description of the spectral projector method for the computa- 
tion of iy{M,n). For a more complete exposition, we refer to the original work of Ref. [30]. 
In this section we assume that we work on an Euclidean lattice, but we will not specify the 
particular form of the lattice Dirac operator. 

If Wm is the orthogonal projector to the subspace of fermion fields spanned by the lowest 
lying eigenmodes of the massive Hermitian Dirac operator D^D with eigenvalues below some 
threshold value M^, the mode number ij{M,fi) can be represented stochastically by: 



\ i=i / 

where rji, . . ., r]j\f are pseudo- fermion fields added to the theory. 

The orthogonal projector IPm can be approximated by a rational function of D^D: 




(2.1) 



X = 1 - 



2M2 



(2.2) 



DW + M2 
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where M* is a mass parameter related to the spectral threshold value M^. The function: 

h{x) = ^{l-xP{x')) (2.3) 

is an approximation to the step function 9{—x) in the range —1 < x < 1, where P{y) is 
in our case the Chebyshev polynomial (of some adjustable degree n) that minimizes the 
deviation: 

5= max |l-^P(y)| (2.4) 

for some e > 0. Computing the approximation to the spectral projector Pjv/ requires solving 
the following equation an appropriate number of times: 

{D^D + M^)i; = 7] (2.5) 

for a given source field rj. Solving this equation is the main computational cost in the 
calculation of the mode number. In particular, the computational cost scales linearly V 
with the volume. 

One can show [30] that the mode number is a renormalization group invariant, i.e.: 

UR{MR,^lR) = u{M,^l), (2.6) 

where the subscript R denotes renormalized quantities. Note that the spectral threshold 
parameter M renormalizes in the same way as the light quark mass (i.e. Mr = Zp^M for 
Wilson twisted mass fermions). 

Finally, we give here the relation between the mode number and the mass-dependent 
renormalized chiral condensate: [30] 



which is defined to match the chiral condensate to leading order of chiral perturbation 
theory. 



3 Lattice setup 

In this section, we will specify the lattice Dirac operator that is used for our work, i.e. the 
Wilson twisted mass Dirac operator. For our computations of the chiral condensate, we 
have used gauge field configurations generated by the European Twisted Mass Collaboration 
(ETMC) with Nf = 2 [6, 32, 33] and A^/ = 2 + 1 + 1 [34-36] dynamical quarks. 
The gauge action is: 

Sg[U] = ^Y.[^^ E ReTr(l-Pi^-^,i)+6i^ReTr(l-P,i^^2^)), (3.1) 

^ As shown in Ref. [30], the ratio M/M* depends on the chosen approximation to the projector. For the 
choice of Ref. [30], which we apply also in our case, M/M* = 0.96334. 
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with j3 = 6/(7o, 5o the bare couphng and P^^^, P^^^ g^j,g ^j^g plaquette and rectangular 
Wilson loops, respectively. For the Nf = 2 case, we use the tree-level Symanzik improved 
action [37], i.e. we set hi = — with the normalization condition = 1 — 8bi. In the case 
of Nf = 2 + 1 + 1, we use the Iwasaki action [38, 39], i.e. hi = —0.331. 

The Wilson twisted mass fermion action for the light, up and down quarks for both 
the Nf = 2 and Nf = 2+1 + 1 cases, is given in the so-called twisted basis by: [40-43] 

= a'^'^Xi{^){Dw + rno^i + il^n5T3)xi{x), (3.2) 

where uiq i and fii denote, respectively, the bare untwisted and twisted light quark masses 
(for shortness, whenever there is no risk of confusion, from now on we will use the symbol 
fi to denote fii). The renormalized light quark mass is given by fin = Zp^fi. The matrix 
acts in flavour space and Xl = iXu, Xd) is a two-component vector in flavour space, related 
to the one in the physical basis by a chiral rotation. The standard massless Wilson-Dirac 
operator D\y reads: 

Dw = ^(7m(V^ + Vp - aV;V^), (3.3) 

where and are the forward and backward covariant derivatives. 
The twisted mass action for the heavy doublet is given by: [42, 44] 

Sh['>P,i^,U] = a^'^Xh{x){Dw + rno^h + ifJ-crlbTi + f^sT3)xh{x), (3.4) 

X 

where mQ fi denotes the bare untwisted heavy quark mass, the bare twisted mass with 
the twist along the ri direction and fis the mass splitting along the rs direction, introduced 
to make the strange and charm quark masses non-degenerate. The mass parameters fLu and 
Us are related to the physical renormalized strange m|j and charm quark masses by 
m'^j^ = Zp^ {fifj + {Z p / Z s) fJ-s) ■ The heavy quark doublet in the twisted basis Xh = iXc Xs) 
is again related to the one in the physical basis by a chiral rotation. 

The twisted mass formulation allows for an automatic 0{a) improvement of physical 
observables, provided the hopping parameter k = (8 + 2amo)~^, where niQ = rriQ^i = mo,/i 
can be chosen, is tuned to maximal twist by setting it to its critical value, at which the 
PC AC quark mass vanishes [41, 45-49]. 

The details of the ensembles considered for this work are presented in Tab. 1 for 
Nf = 2 and Tab. 2 for Nf = 2 + 1 + 1. For both cases, they include 3 lattice spacings 
(from a ~ 0.05 to a ~ 0.085 fm) and up to 5 quark masses at a given lattice spacing. 
The renormalized light quark masses fj,^ are in the range from around 15 to 50 MeV. The 
values of the renormalization constant Zp for different ensembles^ [52-54], used to convert 
bare light quark masses fi and bare spectral threshold parameters M to their renormalized 
values in the MS scheme (at the scale of 2 GeV), are given in Tab. 3, where we also show 
the values of ro/a (in the chiral limit), used to express our results for the condensate as a 
dimensionless product r^T,^^^. Our physical lattice extents L for extracting physical results 

^For Nf = 2 + 1 + 1, the mass-independent renormalization constant Zp is extracted as a chiral limit of 
a dedicated computation with 4 mass-degenerate flavours - see Refs. [55, 56] for details. 
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Ensemble 


/3 


lattice 


afj, 


I2R [MeV] 




K[. 


L [fm 


b30.32 


3.90 


32^ 


X 


64 


0.003 


16 


0. 


160856 


2.7 


b40.16 


3.90 


le^* 


X 


32 


0.004 


21 


0. 


160856 


1.4 


b40.20 


3.90 


20^ 


X 


40 


0.004 


21 


0. 


160856 


1.7 


b40.24 


3.90 


24^ 


X 


48 


0.004 


21 


0. 


160856 


2.0 


b40.32 


3.90 


32^ 


X 


64 


0.004 


21 


0. 


160856 


2.7 


b64.24 


3.90 


243 


X 


48 


0.0064 


34 


0. 


160856 


2.0 


b85.24 


3.90 


243 


X 


48 


0.0085 


45 


0. 


160856 


2.0 


C30.32 


4.05 


32^ 


X 


64 


0.003 


19 


0. 


157010 


2.1 


C60.32 


4.05 


32^ 


X 


64 


0.006 


37 


0. 


157010 


2.1 


C80.32 


4.05 


32^ 


X 


64 


0.008 


49 


0. 


157010 


2.1 


d20.48 


4.20 


48^ 


X 


96 


0.002 


15 


0. 


154073 


2.6 


d65.32 


4.20 


32^ 


X 


64 


0.0065 


47 


0. 


154073 


1.7 



Table 1: Parameters of the Nf = 2 gauge ensembles [6, 32, 33]. We show the inverse bare 
coupling /3, lattice size (L/a)^ x (T/a), bare twisted light quark mass a/i, renormalized 
quark mass /ijj in MeV, critical value of the hopping parameter at which the PCAC mass 
vanishes and physical extent of the lattice L in fm. 



Ensemble 


/3 


lattice 


am 


fii^R [MeV] 






L [fm 


A30.32 


1.90 


32^ 


X 


64 


0.0030 


13 


0. 


163272 


2.8 


A40.20 


1.90 


20^ 


X 


40 


0.0040 


17 


0. 


163270 


1.7 


A40.24 


1.90 


243 


X 


48 


0.0040 


17 


0. 


163270 


2.1 


A40.32 


1.90 


32^ 


X 


64 


0.0040 


17 


0. 


163270 


2.8 


A50.32 


1.90 


32^ 


X 


64 


0.0050 


22 


0. 


163267 


2.8 


A60.24 


1.90 


243 


X 


48 


0.0060 


26 


0. 


163265 


2.1 


A80.24 


1.90 


243 


X 


48 


0.0080 


35 


0. 


163260 


2.1 


B25.32 


1.95 


32^ 


X 


64 


0.0025 


13 


0. 


161240 


2.5 


B35.32 


1.95 


32^ 


X 


64 


0.0035 


18 


0. 


161240 


2.5 


B55.32 


1.95 


32^ 


X 


64 


0.0055 


28 


0. 


161236 


2.5 


B75.32 


1.95 


32^ 


X 


64 


0.0075 


38 


0. 


161232 


2.5 


B85.24 


1.95 


243 


X 


48 


0.0085 


45 


0. 


161231 


1.9 


D15.48 


2.10 


48^ 


X 


96 


0.0015 


9 


0. 


156361 


2.9 


D20.48 


2.10 


48^ 


X 


96 


0.0020 


12 


0. 


156357 


2.9 


D30.48 


2.10 


48^ 


X 


96 


0.0030 


19 


0. 


156355 


2.9 



Table 2: Parameters of the Nf = 2 + 1 + 1 gauge ensembles [34-36]. We show the inverse 
bare coupling /3, lattice size (L/a)'^ x (T/a), bare twisted light quark mass /i^, renormalized 
quark mass fx^R in MeV, critical value of the hopping parameter at which the PCAC mass 
vanishes and physical extent of the lattice L in fm. 

range from 2 to 3 fm (in the temporal direction, we always have T = 2L). To check for 
the size of finite volume effects, we included different lattice sizes for /? = 3.9, a^u = 0.004 
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Nf p g [fm] Zp(MS, 2GeV) rp/a 



2 4.20 0.054 

2+1+1 1.90 0.0863 

2+1+1 1.95 0.0779 

2+1+1 2.10 0.0607 



2 4.05 0.067 



2 3.90 0.085 



0.437(7) 
0.477(6) 
0.501(13) 
0.529(9) 
0.504(5) 
0.514(3) 



5.35(4) 
6.71(4) 
8.36(6) 
5.231(38) 
5.710(41) 
7.538(58) 



Table 3: The values of the lattice spacing a [36, 50], rg/a [34, 50, 51] and the renormal- 
ization constant Zp in the MS scheme at the scale of 2 GeV [52-54], for different values of 
/3 and iV/ = 2 and iV/ = 2 + 1 + 1. 

{Nf = 2) and /3 = 1.9, a^t = 0.004 (iV/ = 2 + 1 + 1). 
4 Results 

In this section, we show our results of the calculation of the chiral condensate. First, we 
illustrate the procedure of extraction of the chiral condensate and discuss the influence of 
the various errors that enter the computation. Then, we analyze finite volume effects in 
our simulations. Finally, we move on to our chiral and continuum extrapolations. 

4.1 Procedure and errors 

We show here how to extract the mass-dependent chiral condensate according to Eq. (2.7), 
illustrating the procedure for ensemble B40.32. Using the spectral projector method, we 
computed the dependence of the mode number on the renormalized spectral threshold 
parameter Mr for 5 values of M^j, from around 2.5 times the renormalized quark mass to 
around 120 MeV. Shortly above the latter value one starts to see deviations from the linear 
regime of ^^{Mji, fiji) vs. M/j (see Appendix A). 

Fig. 1 shows the dependence of the mode number on the renormalized spectral thresh- 
old parameter Mp. The solid line is a linear fit to all 5 points. The slope of this line 
dvR^Mp, fiji) /dMji determines the value of the mass-dependent chiral condensate accord- 
ing to Eq. (2.7). The error of this slope includes two sources: the error of the slope of 
the bare mode number as function of the bare threshold parameter M, du{M,^)/dM 
and the error of Zp needed to convert from bare to renormalized quantities. Although 
di'ji^Mp, fLp) /dMji appears to be constant as a function of Mr within errors, we will take 
its value to be the middle point of the chosen fitting interval, see below for details of the 
fitting intervals considered. Finally, Eq. (2.7) yields, after taking the cubic root^: 



where the first error is the one of the slope di'(M, fi)/dM and the other one comes from 
Zp = 0.437(7) and is dominated by systematic effects - hence, we take it as a systematic 
error of our computation. 



1/3 



0.13372(34)(72), 
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Figure 1: Dependence of the mode number on the renormaUzed spectral threshold pa- 
rameter Mn for ensemble B40.32. The solid line is a linear fit to all 5 points. 



fit 1 


-ange 


lowest oMr 


highest aMn 


xVdof 


roSi/3 


1 


- 3 


0.0225 


0.0375 


0.004 


0.7085(43) 


1 


- 4 


0.0225 


0.0450 


O.OfS 


0.7116(24) 


1 


- 5 


0.0225 


0.0525 


0.588 


0.7154(18) 


2 


- 4 


0.0300 


0.0450 


0.006 


0.7141(46) 


2 


- 5 


0.0300 


0.0525 


0.567 


0.7189(33) 


3 


- 5 


0.0375 


0.0525 


0.549 


0.7235(57) 



Table 4: Values of r^T}^'^ extracted from different fitting ranges. Every fit includes at 
least 3 values of Mr. The fit labeled "1 - 5" is the full fit. We estimate the error from the 
choice of the fitting range by comparing the value from the full fit with the ones from fits 
"1 - 4" and "2 - 5". The error given is statistical only. 

The value of aS^/^ can be further converted to a dimensionless product rgS^/^ (which 
will be the final result of this paper, after taking the chiral and continuum limits) or to a 
physical value in MeV. For the former, we use the value in the chiral limit r^/a = 5.35(4). 
Since the error of this value is again mostly systematic, we quote it as another systematic 
error of tqT}/^: 

= 0.7154(18)(39)(53), 
^The values of E that we give (also in our plots) are always for the renormalized condensate. 
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where the two errors are as above and the third one comes from ro/a. For a conversion to 
MeV, one needs to choose a value of the lattice spacing in physical units. There are several 
such estimates for ETMC 2-flavour ensembles, giving for /3 = 3.9 values including 0.079 fm 
[6], 0.085 fm [50] and 0.089 fm [57]. Taking this spread into account, the relative error on the 
lattice spacing is around 7%, which leads to a similar relative uncertainty in the value of the 
chiral condensate S^/^ in physical units, which amounts to about 20 MeV. This is roughly 
an order of magnitude more than other errors entering our computation. Even being less 
conservative and using for the error the value quoted in Ref. [50] - 0.085(2)stat(l)syst ^ 
the error that it yields is still of the order of 10 MeV. Therefore, we decided to give our 
final results as the dimensionless product roS^/'^ and we chose not to quote any value in 
MeV for it until a significantly improved determination of the lattice spacing is available. 

Another source of the error is the choice of the fitting range and hence the value of 
M/j that enters the square root in Eq. (2.7). Of course, physical results should not depend 
on this choice, provided that the whole fitting interval lies in the linear regime of the 
mode number vs. Mr dependence. Hence, varying the fitting range serves two purposes: 
establishing whether non-linear effects are already present and checking that the choice of 
Mr in Eq. (2.7) does not influence the final result. The values of r^T}/'^ resulting from 
different fitting ranges in Mr are shown in Tab. 4. For all fits, x^/d-O-f. is below 1. The 
compatibility of all results and the good values of x'^/d.o.f. imply that for this ensemble 
the choice of the fitting range and Mr does not affect the results in a substantial way^. To 
quantify this error, we considered the 4-point fits with the lowest or highest value of Mr 
excluded. Excluding the first or last point leads to a similar change of the result and hence 
we took the larger of the two as our conservative error from the choice of the fitting range. 
Note, however, that Tab. 4 implies a systematic tendency towards increasing of S when 
the fitting range moves towards higher values of Mr. This indicates an onset of non-linear 
behaviour for values of Mr only slightly above the ones we considered. 

Finally, our estimate of r^T}/^ for ensemble B40.32, including all sources of error, is: 

= 0.7154(18)stat(38)fitrangc(39)zp(53),„/„. 

We note that the total error is dominated by systematic errors. This means that increasing 
statistics would not essentially change our total error. It should be considered an important 
advantage of the method of spectral projectors that rather moderate statistics (in our case 
around 230 independent gauge field configurations for this ensemble) leads to a practically 
negligible statistical error. 

4.2 Finite volume effects 

One of the main sources of systematic effects in Lattice QCD simulations are finite volume 
effects (FVE). In Ref. [30], theoretical arguments were provided that FVE should be small 
for the chiral condensate computed from the mode number - with exponentially small 
difference between the finite volume and infinite volume results of 0(exp(— MaL/2)), where 

*The ensemble B40.32 is somewhat special in this aspect. As we show below, in general, the fitting range 
uncertainty is the most important source of error in our analysis. 
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aMR=0.0225 
-aMp=0.0300 I — - 
aMR=0.0375 
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aMR=0.0525 I — - 
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20 24 
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32 



Figure 2: The main plots show the vohime dependence of the mode number density v/V 
for different values of the renormalized spectral threshold Mr. The horizontal bands show 
the result at the largest volume. The insets show the volume dependence of the chiral 
condensate rgS^/^ (the error of each point includes the statistical error and the systematic 
one originating from the choice of the fitting interval). (Left) Nj = 2, /3 = 3.9, a/x = 0.004, 
(right) Nf = 2 + 1 + 1, /3 = 1.9, a// = 0.004. 



= 2A5]/F^, A = — /i^ and F is the pion decay constant in the chiral limit. Since 

in practice the mass-dependent chiral condensate is extracted at A ^ ^, the mass M\ is 
much higher than the pion mass, which typically governs FVE. Hence, one expects that for 
the computation of the chiral condensate from the mode number, FVE will be rather small. 
FVE for the mode number itself were computed in SU(2) chiral perturbation theory [23]. 
The resulting formula leads to a prediction that FVE from lattices with L > 2 f m should be 
smah, 0{< 1%) for Mr ^ 0(60 - 120) MeV and renormalized quark masses of 0(10 - 20) 
MeV (with larger FVE at smaller Mr)- Indeed, in practice, it was shown in Ref. [30] that 
for L > 2 fm the results deviate from their infinite volume values by less than 1%. 

To show that it is also the case in our setup, we performed the computation of the 
mode number and the chiral condensate for: 

• Nf = 2 : 4 different volumes at fixed /3 = 3.9, a/i = 0.004, lattice extents: L/a = 
16, 20, 24 and 32, with corresponding physical values of 1.4, 1.7, 2.0 and 2.7 fm, 
respectively, 

• Nf = 2 + 1 + 1 : 3 different volumes at fixed /3 = 1.9, a/i = 0.004, lattice extents: 
L/a = 20, 24 and 32, with corresponding physical values of 1.7, 2.1 and 2.8 fm, 
respectively. 

In Fig. 2, we show the volume dependence of the mode number density u/V for 4-5 
different values of the renormalized spectral threshold Mr. The mode number density can 
be computed very precisely. It hence provides a strong test of finite size effects. 

The left plot shows our data for Nf = 2. The results for L/a = 20 and especially 
L/a = 16 are systematically lower than L/a = 32, signaling large FVE. However, the mode 
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number density for L/a = 24 is compatible with the one for L/a = 32 for 3 intermediate 
values of Mr, while it differs by 2-3cr for the lowest and highest M^j, thus changing the slope 
of the mode number vs. dependence and the extracted chiral condensate (see the inset 
of Fig. 2 (left)). This change of slope is statistically significant, but it is still a relatively 
small, 1-1.5% effect. Taking into account the uncertainty from the choice of the fitting 
range, the final results for the chiral condensate are compatible for all cases - including the 
ones for small volumes, indicating that even if the mode number density goes systematically 
down, the slope of the whole z^(Mr,^/j) dependence is less affected. We have also tried a 
description of FVE in the framework of the formula derived in Ref . [23] . We conclude that 
it provides a reasonable agreement with actual lattice data for L/a > 24, while FVE for 
smaller volumes are somewhat underestimated (by a factor of 0(2) at L/a = 16, compared 
to the actually observed FVE). 

In the right plot, we show analogous data for Nf = 2 + 1 + 1. Similarly, we observe 
significant finite size effects in the mode number density (and also in the chiral condensate) 
for L/a = 20, while L/a = 24 and L/a = 32 are always compatible. 

This allows us to conclude that finite size effects are indeed very small when one reaches 
a linear lattice extent of around 2 fm {L/a = 24 at both f3 = 3.9 {Nf = 2) and /3 = 1.9 
{Nf = 2 + 1 + 1)). Therefore, we used in our analysis all available ETMC ensembles with 
a linear lattice extent of at least 2 fm^. 

4.3 Chiral and Continuum Limit — Nf = 2 

We now show our results for the 2-flavour case. For each value of (3, we have 2-4 sea quark 
masses, according to Tab. 1. For each ensemble, we perform computations of the mode 
number at 5 values of the renormalized spectral threshold Mr, from around 50 to 120 MeV. 
We follow the procedure outlined in Sec. 4.1, i.e. we extract the mass-dependent condensate 
from the slope ^{{{Mr, IJ-r) as function of Mr for each ensemble. 

The results for rgS^/^ for all considered ensembles are gathered in Tab. 5. These results 
are then used to extrapolate to the chiral limit for each value of f3. Following the argument^ 
outlined in Ref. [30], we perform a linear extrapolation of a^T, vs. a/j, to the chiral limit. 
Our extrapolations for all three values of /3 are shown in the main plot of Fig. 3 (we plot 
TqS vs. rQfj,f( to allow comparisons between different values of /3). The plotted errors are 
only statistical, since in extrapolations at fixed /3, the relative errors from Zp and rg/a are 
the same for all quark mass values (we use chirally extrapolated values of Zp and ro/a) - 
we give them in Tab. 5. To estimate the systematic error originating from the choice of 
the fitting range in vp{Mp, ^p) vs. Mp fits, we repeated all chiral extrapolations for two 
tailored fitting ranges - excluding the first value of Mp (to account for effects of coming 

^The exception to this rule is ensemble d65.32 with L ~ 1.7 fm. We decided to use this ensemble, 
because without it there would only be one quark mass at /3 = 4.2 and no chiral extrapolation could be 
performed. 

''The next-to-leading order chiral perturbation theory expansion of the mass-dependent condensate sug- 
gests that the latter is equal to the chiral condensate in the chiral limit up to terms linear in and higher 
order effects. In particular, there are no corrections proportional to In^/j. 
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Table 5: Results for rgS^/'^ for all considered Nf = 2 ensembles. The given errors are: 
statistical, resulting from Zp, resulting from ro/a, respectively. We also show results in the 
chiral limit, where we also give the systematic error from the choice of the fitting interval 
(4th error). See text for more details. 



too close to the renormalized quark mass) or the last value thereof (to account for possible 
deviations from the linear behaviour for too high values of M/j). 

The chiral limit values, with all sources of error, are also shown in Tab. 5. In general, 
the total error originates in practice only from the choice of the fitting range and the latter 
increases when approaching the continuum limit. The reasons for this behaviour include the 
fact that the number of quark masses that we use decreases for smaller lattice spacings and 
at /3 = 4.2 the slope of the quark mass dependence of the chiral condensate is apparently 
larger than at coarser lattice spacings^, making the final chiral limit value more susceptible 
to changes in the fitting interval. 

Finally, we can use the chirally extrapolated values of the condensate to perform an 
extrapolation to the continuum limit. The chiral condensate is an C'(a)-improved quantity, 
i.e. it is even under the TZ^ parity transformation: ^p — ?■ i'y^T^ip, 'ip — ?• i'il)^^T^ [41]. This 
can be shown by introducing the spectral sums [30]: ak{n,mq) = (Tr ^{DmDm + //^)~'^|), 
where k > 3 for reasons explained in the given reference. The spectral sums are related to 
the mode number [30] and the improvement of the spectral sums implies the improvement 
of the chiral condensate. Representing the spectral sums as a density chain correlation 
function (for k = 3): 

as{fi,m) = -a24 ^ (P+(xi)P2-3(x2)P34(^3)^75(^4)P5'6(^5)^6"l(^6)), (4.1) 
XI. ..xe 

it is straightforward to show that the object on the right-hand side is even under TZ^ 

^Note that this slope may be affected by the smaller volume of ensemble d65.32 - hence, it may be a 
residual FVE and not an indication of the dependence of the slope on the lattice spacing. In such case, our 
fitting range error at /? = 4.2 implicitly reflects this FVE. 
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0.02 0.04 0.06 0.08 0.1 0.12 

Figure 3: Main plot: chiral extrapolations of the chiral condensate TqE for Nf = 2 
ensembles and /3 = 3.9, 4.05, 4.2. The lines are extrapolations to the chiral limit, linear in 
the quark mass. The values in the chiral limit for /3 = 3.9 and 4.2 are slightly shifted for 
better presentation. The errors are statistical only. Inset: continuum extrapolation of the 
chirally extrapolated chiral condensate roS^/^ vs. (a/ro)^. The errors include: statistical 
errors, errors from Zp and errors from rg/a. 

transformation, which changes the sign of the charge to the individual pseudoscalar densities 
and the overall sign. However, since we always have an even number of densities, the latter 
remains positive, yielding C'(a)-improvement of the density chain and hence the chiral 
condensate. For a complete proof of C'(a)-improvement, one also needs to consider contact 
terms arising from Eq. (4.1), i.e. terms in the sum with Xi = xj for some i ^ j, which 
can potentially give rise to terms linear in the lattice spacing. We are currently working on 
completing this proof. Nevertheless, we observe relatively small cut-off effects in the chiral 
condensate and hence our conclusions will remain essentially unchanged even if such linear 
terms arise, i.e. the extracted continuum limit values assuming 0{a) scaling are consistent 
(although with roughly 50% larger errors) with the ones assuming 0{a'^) leading cut-off 
effects. 

The continuum limit extrapolation is performed linearly in a^, using results at three 
lattice spacings, with fixed fitting range of the mode number vs. spectral threshold depen- 
dence, corresponding to Mr ~ 90 MeV (entering the square root in Eq. (2.7)) for all values 
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of (3. As an error, we use the statistical error, combined in quadrature with the error of Zp 
and ro/a. We do not observe significant cut-off effects. The final value in the continuum 
limit is 0.689(16). To account for the fitting range error, we perform the full analysis for 
tailored fitting ranges, excluding the first or last value of Mp for each ensemble. This corre- 
sponds to a shift in Mp to approx. 80 or 100 MeV, respectively. While the extracted value 
of the condensate in the chiral limit should not depend on the fitting range, in practice the 
results for different fitting ranges differ, which is due to using only 4-5 values in the fits to 
extract the slope of ^^{Mii, ijlr). The fits from tailored fitting ranges yield 0.678(18) and 
0.718(20), respectively. To be conservative, as our systematic error from the fitting range 
we choose the larger difference of the two with respect to the central value 0.689(16). This 
finally gives: 

roS;^'^2 = 0-689(16) (29), 

where the first error is the combined statistical error, the error of Zp and of ro/a, while the 
second error originates from the choice of the fitting range. 

4.4 Chiral and Continuum Limit — Nj = 2 -|- 1 -|- 1 

In this subsection, we present results for the Nj = 2 -|- 1 -|- 1 case. By comparing to the 
results of the 2-flavour case, we can investigate the role of the dynamical strange and charm 
quarks. 

We proceed as in the previous section. For each value of f3, we have 3-5 sea quark masses, 
according to Tab. 2. We compute the mode number at 4 values of the renormalized spectral 
threshold Mr, from around 50 to 110 MeV, and extract the mass-dependent condensate 
from the slope of the iyji{Mp, ^p) vs. Mp dependence for each ensemble. We have also 
computed the mode number for a fifth value of Mp ~ 130 MeV. However, given the very 
good statistical precision of the spectral projectors method of evaluating the mode number, 
we observe significant deviations from the linear dependence of i^p{Mp, fip) on Mp when 
this fifth value of Mp is included. Because of this, we decided not to include the results at 
Mp w 130 MeV. 

In Tab. 6, we show all our results for rgS^/'^ in the 2+1 + 1-flavour case. We also include 
the results of a linear extrapolation to the chiral limit for each value of /3, shown in the 
main plot of Fig. 4. As before, we plot only statistical errors, since all extrapolations are 
performed at fixed P and the errors from Zp and ro/a are the same for all quark masses 
(we use chirally extrapolated values of Zp and ro/a) - given in Tab. 6. Contrary to the 
Nf = 2 case, the slope of the dependence of the mass-dependent condensate on the light 
quark mass slightly decreases for increasing /3 (the change of slope is statistically significant 
when going from /? = 1.9 to /? = 2.1). This has the effect of lowering the systematic error 
related to the choice of the fitting range for decreasing lattice spacing^, which is for all (3 
the dominating source of error (although for /? = 2.1 other errors become comparable). 

The chirally extrapolated values at three lattice spacings are then used to perform an 
extrapolation to the continuum limit, which is again compatible with 0{a^) cut-off effects. 

^Moreover, we always have at least 3 sea quark masses for each /3 in the Nf = 2 + 1 + 1 case, compared 
to only 2 masses at /3 = 4.2 {Nf = 2)). 
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Table 6: Results for roT}/'^ for all considered Nf = 2 + 1 + 1 ensembles. The given errors 
are: statistical, resulting from Zp, resulting from ro/a, respectively. We also show results 
in the chiral limit, where we also give the systematic error from the choice of the fitting 
interval (4th error). See text for more details. 



To estimate the fitting range uncertainty, we again perform 3 separate continuum limit 
extrapolations, using different fitting ranges and different values of M/j, corresponding to 
approx. 80, 90 and 100 MeV. The values of roY^I^ in the continuum limit are, respectively, 
0.668(24), 0.680(20) and 0.659(27). As our central value we take the result from the fuU 
fitting range: 

roS;^^L2+i+i = 0.680(20)(21), 

with the larger of the differences with respect to values from tailored fitting intervals as the 

fitting range systematic error. This can be compared to the 2-flavour result which amounts 
1 /3 

to ro5]j^^,^2 = 0.689(16) (29) and both results are compatible within errors. 
5 Conclusions 

In this paper, we presented our results on the chiral condensate in QCD with Nf = 2 and 
Nf =2 + 1 + 1 fiavours of dynamical Wilson twisted mass quarks at maximal twist. 
Our final results are: 



^oS;^'=2 = 0-689(16) (29), 

.1/3 



roS;;;^2+i+i = 0-680(20) (21), 

which indicates that at the current level of precision, we cannot discriminate the infiuence of 
the dynamical strange and charm quarks on the value of the light quark chiral condensate. 
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Figure 4: Main plot: chiral extrapolations of the chiral condensate TqE for Nf = 2 + 1 + 1 
ensembles and (3 = 1.9, 1.95, 2.1. The lines are extrapolations to the chiral limit, linear 
in the quark mass. The errors are statistical only. Inset: continuum extrapolation of the 
chirally extrapolated chiral condensate rgS^/^ vs. (a/rof. The errors include: statistical 
errors, errors from Zp and errors from rg/a. 



The main source of the error of our results is the systematic error related to the choice 
of the fitting range in the dependence of the renormalized mode number on the renormalized 
spectral threshold. The second most important source of the error is either the statistical 
error or the error related to the uncertainty in the values of ro/a (which are inputs of our 
analysis). The error from Zp (also an input of our analysis), used to renormalize the quark 
masses and the spectral threshold parameter, is usually the smallest. However, we want to 
emphasize that in all cases the fitting range error is the largest one and in most cases it is 
larger by a factor of 2-4 than any other error. The rather small statistical errors that we 
obtain indicate that increasing statistics would not make our total error significantly smaller. 
This implies that a way to improve the total error would be to increase the number of values 
of the spectral threshold Mp at which one computes the mode number. This would allow 
to identify more precisely the linear region of the mode number vs. Mp dependence (see 
discussion in the Appendix) - on the one hand sufficiently far away from the renormalized 
quark mass, on the other hand low enough such that there are no deviations from the linear 
behaviour at the upper end of the fitting range (observed already at Mp ~ 130 MeV in the 
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Result 



method 



N 



f 



fermions 



this work 
this work 



spectral proj. 2 twisted mass 0.689(16)(29) 

spectral proj. 2+1 + 1 twisted mass 0.680(20)(21) 



RBC-UKQCD [3] 
MILC [4] 
MILC [5] 
!. Borsanyi et al 
ETMC [6] 
ETMC [26] 
HPQCD [27] 



chiral fits 2+1 

chiral fits 2+1 

chiral fits 2+1 

[9] chiral fits 2+1 

chiral fits 2 

quark propagator 2 

quark propagator 2+1 + 



1 



domain wall 
staggered 
staggered 
staggered 
twisted mass 
twisted mass 
staggered 



0.632(15)(12) 
0.654(14)(18) 
0.653(18)(11) 

0.662(5)(20) 
0.575(14)(52) 
0.676(89)(14) 

0.673(5)(11) 



Table 7: Comparison of our results for roS^/'^ with large-volume continuum limit results 
found in literature. In the given references, the values of E are given in MeV. To convert 
to the dimensionless product r^T,^/'^, we combine them with results for rg. The first error 
given is always from the computation of the value in MeV (when there are several errors 
given, we combine them in quadrature) and the second one from conversion using physical 
value of ro. For RBC-UKQCD, = 256(6) MeV, tq = 0.487(9) fm [3]. For MILC [4], 
SVs = 278(6) MeV, n = 0.318(7) fm [4], ro/ri=1.46(l)(2) [58]. For MILC [5], S^s = 
281.5(7.9) MeV, n = 0.3133(23) fm [5], ro/ri=1.46(l)(2) [58]. For S. Borsanyi et al. [9], 
SV3 = 272.3(1.2)(1.4) MeV, ro = 0.48(1)(1) fm [59]. For ETMC [6], S^/^ ^ 269.9(6.5) 
MeV, ro = 0.420(14) fm. However, newer analyses indicate a higher value of tq ~ 0.45 fm 
[50]. To take this into account, we added the spread of the new and old value as a systematic 
error and used ro = 0.420(38) fm to calculate roS^/^ f^-^j^ -^^ ^^y^ p^^. g^MC [26], 
= 299(26) (29) MeV, ro = 0.446(9) fm [6]. For HPQCD [27], S^s = 283(2) MeV, 
ri = 0.3209(26) fm [60], ro/ri=1.46(l)(2) [58]. 



Nf = 2 + 1 + 1 case). 

In order to place our values of the chiral condensate in context of results of other col- 
laborations, we attempt in Tab. 7 a comparison. Given the large amount of approaches 
to compute the chiral condensate, as mentioned in the introduction, we make a selection 
by only considering results that are given in the literature as continuum limit values from 
large volume simulations. Note that the other available results that are listed in Tab. 7 
are obtained in a different way than the spectral projector method. They mostly use chiral 
perturbation theory fits to the quark mass dependence of light pseudoscalar meson observ- 
ables or determine the chiral condensate from the quark propagator. Not discussing the 
advantages or disadvantages of the various fermion discretizations used, we see in Tab. 7 an 
overall agreement for the dimensionless quantity roS^/^, which is reassuring and establishes, 
in our opinion, the spectral projector method as a valuable alternative to determine the 
chiral condensate. On the other hand, when looking at the chiral condensate in physical 
units, see the caption of Tab. 7, a spread of results is obtained. Thus, it seems that the scale 
setting from the different lattice calculations introduces a systematic effect and it would be 
desirable to clarify this uncertainty in future more precise calculations. 
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A Testing the method 

To test our implementation of the method in the tmLQCD package [61], we compared the 
spectral projectors results for the mode number with ones from explicit computation of 
150 lowest eigenvalues of D^D for each gauge field configuration, using ensemble B85.24. 
The results of this comparison are shown in the right plot of Fig. 5, where the 4 points 
correspond to the stochastically evaluated mode number using spectral projectors, while 
the continuous line (which has an error roughly of the order of the width of the line) is the 
result from explicitly computing the eigenvalues. We observe very good agreement between 
the two methods. 

Moreover, we used the results of explicit computation of eigenvalues to estimate the 
region of renormalized spectral threshold of Mr where we observe linear dependence be- 
tween the renormalized mode number and the threshold value (left plot of Fig. 5). The 
onset of non-linear behaviour corresponds to approx. 130-150 MeV. On the other end of the 
spectrum, one clearly observes effects of Mr close to the renormalized quark mass up to 
around 10-20 MeV above the latter^. This allows us to identify the linear region to extend 
between around 60 and 120 MeV, to allow for some safety margin. Therefore, we decided 
to choose our values of Mr for the computation of the chiral condensate roughly in this 
interval. We remark that values in this range were used in Ref. [30]. 

Some parameters employed in the method of spectral projectors need to be tuned to 
obtain a compromise between the accuracy of results and computational cost. 

First of all, as mentioned in Ref. [30], the precision of the inverter can be chosen to be 
relatively sloppy without reducing the accuracy. In order to identify the optimal precision, 
which does not affect the correctness of the result, but still decreases the computational 
time, we computed the mode number for several values of the relative precision of the 
inverter. The results (for ensemble b40.16) are shown in Fig. 6, which shows that even 

^One expects that the behaviour close to the renormahzed quark mass is modified in an important way 
by lattice artefacts [23]. 
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Figure 5: The mode number as a function of M/j for the ensemble B85.24. The solid line 
in both plots is the result of explicit computation of 150 eigenvalues for each gauge field 
configuration. The right plot is the zoom of the left part, with added results of spectral 
projectors calculation of the mode number for 4 values of Mjj. 



Nf = 2, /3 = 3.9,16^ X 32 



a; 
-a 




Figure 6: The mode number vs. the number of averaged stochastic sources using different 
values of the inverter relative precision for a 16"^ x 32 lattice at /3 = 3.9, afj, = 0.004 
(ensemble b40.16). 



precision of 10~^ gives reasonable result. However, we decided to be conservative and we 
chose a value of 10~® for the relative precision of the inverter. 

We also checked the dependence of the mode number on the number of stochastic 
sources used for each gauge field configuration, shown again in Fig. 6. We observe that 
all results are compatible within error, which matches the suggestion of Ref. [30] that one 
stochastic source should be enough. However, we observe that adding a second source 
might help to considerably reduce the statistical error, which may be important for shorter 
Monte Carlo runs, when the available number of independent gauge field configurations is 
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rather small. Adding further sources does not change the error considerably, because of 
correlations between results obtained from the same gauge field configuration. 
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